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Abstract 

Free energy calculations of two crystalline phases of the molecular compound Sg were performed via 
molecular dynamics simulations of these crystals. The elemental sulphur Sg molecule model used in our MD 
calculations consists of a semi-flexible closed chain, with fixed bond lengths and intra-molecular potentials 
for its bending and torsional angles. The intermolecular potential is of the atom-atom Lennard- Jones type. 
Two free energy calculation methods were implemented: the accurate thermodynamic integration method 
proposed by Frenkel and Ladd ijj] and an estimation that takes into account the contribution of the zero 
point energy and the entropy of the crystalline vibrational modes to the free energy of the crystal. The 
last estimation has the enormous advantage of being easily obtained from a single MD simulation. Here 
we compare both free energy calculation methods and analyze the reliability of the fast estimation via 
the vibrational density of states obtained from constrained MD simulations. New results on a- and a'-Sg 
crystals are discussed. 
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I. INTRODUCTION 



In recent papers we studied the crystalline phases of the molecular compound Sg, the most 
abundant natural form of elemental sulphur around ambient pressure and temperature (STP). Three 
crystalline allotropes of this compound are known: a-Sg phase is orthorhombic with 16 molecules 
in the non-primitive unit cell lO, 3]- If several samples of this crystal are slowly heated, some of 
them will show a phase transition to jS-Sg at 369 K, which in turn melts at 393 K but most of 



them show a metastable melting point at 385.5 K M,\M- Monoclinic p-Sg is usually obtained from 
the melt lO] and has 6 molecules per primitive unit cell, two of them orientationally disordered 
above 198 K [4]. Monoclinic y-Sg has four molecules in a pseudo-hexagonal close packed unit 



cell, with a density 5.8% higher than a-5g at STP 

Using a very simple molecular model consisting of a cyclic semi-flexible chain of 8 atoms, 
with constant S-S bond lengths, it is possible to reproduce many features of the complex phase di- 
agram of the 5g molecule laBl- The intra-molecular potential model includes a harmonic bending 
potential for S-S-S angles and a double well for the torsional angles jsl, the intermolecular poten- 
tial used was a simple Lennard-Jones non-bonded atom-atom interaction, the details and potential 
parameters are given in the Section Hill Using this simple interaction model we could reproduce, 
via a series of classical constant temperature- constant pressure simulations flQ], the following 
experimental facts: the crystalline structure, the configurational energy and the lattice, bending 
and torsional dynamics (as given by the calculated density of vibrational modes) of a-, P- and y-Sg 
for T>200 K; the orientational dynamical order - disorder phase transition of P-5g and, finally, 
the solid-liquid phase transition of a cubic disordered sample was calculated near the experimental 
value |E,0]. 

Nevertheless, we found a fact that cannot be reproduced by this simple molecular model I^Q]: 
when the temperature of a a-5g MD sample of 288 molecules is lowered below 200 K, our or- 
thorhombic a-Sg sample shows a structural phase transition to a monoclinic phase, with a molec- 
ular array similar to that of a-5g and that we called a'-5g This has not been experimentally 
observed. This spontaneous change was most probably due to the large fluctuations associated to a 
relatively small sample. A larger sample of 512 molecules didn't show this distortion, but its con- 
figurational energy and volume per molecule remained, nevertheless, higher than those calculated 
for a'-Sg iQ]. Fig. [l] shows both structures and the configurational energy vs. temperature, that 
of <x-S% was calculated by slowly decreasing T in steps of 25 to 50 K. The same study for a'-5g 



is performed by increasing T, after a long stabilization and annealing to obtain a totally ordered 
structure at the starting point of 50 K, this is the reason for the different Uconf values, in fig. [H at 
this point. Experimental measurements of the a-Sg structure at 300 and 100 K clearly disregard 
a structural phase transition, unless a metastable state has been experimentally measured 

m. 

and this fact shows the limits of applicability of the simple molecular model. 

In refs. we analyzed if an anisotropic non- bonded atom-atom intermolecular potential 

n 

model was able to reproduce all the experimental facts. Using the program GAMESS |12], we 
performed ab initio calculations of the electronic density distribution of the molecule, and 
the measured anisotropy of this distribution around each atom due to the atomic lone pairs. The 
calculated electronic density reproduces the experimental crystalline measurements of Coppens 
et. al. with a lone pair center at 0.7A of the S location; nevertheless, at the nearest crystalline 
distances found between non-bonded atoms, the atomic anisotropy is low. The MD simulations 
performed with this anisotropic potential did not irnprove the previous result (shown in Fig. [U), 
unless a quite unrealistic atomic anisotropy is used lllifl^ . 

It has to be pointed out that the search of a reliable molecular model for elemental sulphur 
molecules is extremelly useful due to the practical impossibility of performing quantum mechani- 
cal simulations of the complex phase diagram of these molecular crystals, with a large number of 
atoms in the primitive cells. 

In this paper we decided to implement a reliable calculation of the free energy of a-Sg and a'-Sg 
crystals in order to review the conclusions obtained with the first molecular model I^E], based on 
an estimation of the free energy that takes into account the contribution of the vibrational modes 
to the zero point energy and entropy of the sample flQ]. Our present calculations are useful to 
check our previous free energy estimations and, at the same time, for the comparison of both free 
energy calculation methods. The usefulness of this comparison is due to the fast way in which the 
free energy can be estimated, within the quasi-harmonic approximation, using the data of a single 
MD simulation. 

In the following sections we give the details of the implemented free energy calculations, the 
inter- and intra-molecular potential model used, the performed MD simulations, the obtained re- 
sults and our conclusions. 
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Figure 1: Configurational energy of a and a'-S^ as obtained in ref. Q]. The structural transition to a'-S^ is 
observed spontaneously below 200K. 



II. FREE ENERGY CALCULATIONS 



In order to obtain reliable free energy values we implemented the thermodynamic integration 
method of Frenkel and Ladd [Q], including the two correction terms given in refs. Iislfl^ . From 
the results of a single MD simulation, we also estimated the free energy in the quasi-harmonic 
approximation, that is usually employed in lattice dynamics calculations IitI . In Section HVl we 
give the details of our performed MD simulations. 



A. Thermodynamic integration method 

This method proposes a thermodynamic integration along a reversible path between the system 
we are interested in and another one for which the free energy can be analytically calculated 
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yj, US UM- For crystalline samples the reference system is an Einstein crystal, for which the 
internal energy is calculated as: 

N 

UEinstein = U {{r^Q}) + £ a(r, - ro)2 , 

(=1 

where {iq'} is the set of coordinates of the minimum energy configuration (the equilibrium atomic 



coordinates) and U ({r^}) its potential energy, the second term is the harmonic potential energy 
of the N oscillators. The force constant a is taken equal for all degrees of freedom and its value, 
a = 18.0844 kJ/mol, is chosen so as to obtain a mean square displacement equal to the average one 
in the real crystal: {^r^) Einstein ~ i^^^)reai- potential energy difference between the real and the 
reference crystals can be calculated along a reversible artificial pathway that links both systems. 
At any intermediate point of the path, given by the coordinate X {0 <X< 1), the local value of U 
is 

6^(X) = t/(K}) + (l-X)[t/({i^})-f/(K})]+X£a(r,-ro)2. (1) 

(=1 

So, for X = the potential energy of the real crystal is recovered and, for ?i=l the potential 
energy is that of the Einstein crystal: 

0(1 = 0) = Usist 

and 

U{X= 1) = U Einstein ■ 

The free energy difference between the real and the reference crystal is given by 



Q dX Jo dX 



and deriving eq. [T] 



F = FEinstein + J^_^dl(^f^^ a(r, - ro)^ - [[/(r^) - ^ . (2) 



^'=1 / X 

The free energy of the Einstein crystal, with a fixed center of mass for the MD sample, is 



^^U(rS)^'-^H^)-'-^-'^+m(A), (3) 
^ ^' 2^(3 ^ap^ 2^ ^ ^ ^ 

where A = t '^^ Broglie thermal wavelength , P = is the number of atoms 

and Vq is the equilibrium volume of the sample. 

If N,noi is the number of molecules of the system, the final expression, including a correction 
term due to the center of mass constraint of the sample UM and a second correction term 

due to its finite size is: 



The accuracy of this calculation is detemiined by a correct evaluation of the last term in eq. HI 
where long MD simulations are needed to measure the integrand for each value of A,. 

For all values of X, the center of mass of the sample is held fixed at the origin because the 
'springs' of the reference Einstein crystal are tied to the equilibrium atomic positions. This is 
important to avoid divergences of the integrand in eq. lUwhen X is close to I1II16I. 

The third term in eqs. |3] and |4l includes the change of the sign mentioned in ref. Iisi . with 
respect to the original work of Frenkel and Lad d lUj] , due to a new consistent calculation of the 



partition function of the constrained CM system Ill5ll . This term is negligible in a large system but 
must be taken into account to compute absolute free energies of a typical MD sample. 

B, Quasi-harmonic approximation 

Following Born and Huang IitI . the free energy of a harmonic crystal with Nd degrees of 
freedom ( i. e. independent harmonic oscillators) can be calculated as: 

F = £f; = C/,„„/+-£/zv,+fc6r£ln(l-e (5) 

^ (=1 i=l 

where Uconf is the potential energy and the V/ are the vibrational frequencies of a system with 
Nd degrees of freedom. 

We have implemented this method, as in ref. |Q,Q], in two steps. First we calculate the atomic 
velocities self correlation function C{t), using the MD data stored in a free trajectory in the phase 
space of the sample: 

^ mm) 

(v(O).v(O))' 

where <> imply averages over atoms and over different initial times. 

In a second step we calculate the vibrational density of states D{y) via a Fourier transform of 
C{t), and replace the sums over frequencies in eq. |5]by an integral on frequencies, weighted by 
the calculated density of vibrational modes. We have to take into account that the normalization 
factor for D(v) is such that / D{v) dv — Nd- In our case is Nd = l6Nmoi, because we have 8 bond 
constraints per molecule. 



We can then calculate the internal energy, the entropy and free energy as 111 



If f _ hv 

F = Uconf + 2 y dvD{v) hv + ksT J JvD(v) ln(l - e ^) , (6) 



hv 

—^ D{y)dv 



hv 



If f e '^B^ h\ 
U^F-TS = Uconf + / dvD{v)hv+ / -^—D{v)dv. 

The approximations involved in our quasi-harmonic calculations, with the results of a single 
MD simulation, are the following: 

• Uconf is calculated by performing averages, over time and molecules, of the potential energy 
per molecule calculated in the MD run, Uconf =< U >, it is not the value of the minimum 
potential energy Uq of the Einstein crystal in equation|3l 

• The calculated vibrational density of states D{\) of the MD crystalline sample, due to its 
finite size ( a few number of primitive unit cells are included in it) does not give an accurate 
measurement of the density given by the real dispersion curves, the frequencies are only 
measured in a few points of the reciprocal space. 

• The anharmonic frequencies of the sample can be accurately measured with a MD simula- 
tion, but eq. |6lis called the 'quasi-harmonic' approximation because eq. |5]is exact only for 
harmonic potentials. 



III. INTRA- AND INTERMOLECULAR POTENTIAL MODEL 

The potential model is that of refs. Q]. The flexible molecular model includes all low 
frequency molecular distortions that mix with lattice modes and can therefore be of relevance in 
a possible structural phase transition. The S-S bond distances are kept constant {dss = 2.0601 
A) because the stretching modes are well above in energies (v > 400 cm^^) than the rest of the 
vibrational modes (v < 250cm^^). 

We must stress here that bond constraints can have nonnegligible contributions to the free 
energy calculations. Such contributions have been studied elsewhere , with different methods and 



techniques 1 19 , 



2C 



2 ID . The main conclusions of those explicit calculations, using both MD I19I1 
and MCTI Il20l] . are that the influence of bond constraints is actually negligible when the bond 
length is not changed from the initial to the final state. This is the case in our Sg phases, in which 
the primitive cell will change and the molecules will distort from one phase to the other but the 
molecules bond length is held fixed within the same value for every calculation presented in this 
work. If changes in bond length between initial and final states of the free energy calculation were 
involved, an explicit constraint contribution to the free energy must be considered as explained in 
refs. I^-Q]- For the sake of completeness, a simple estimation of the change in the probability 
density between the standard NVT and the 'bond-constrained' ensemble of is given in the 
Appendix. 

The bending intramolecular potential for S-S-S angles is harmonic 
with a force constant of Cf!=25725 fce/rad^ and |3o= 108 deg. The intramolecular potential for 

A 

torsion angles is a double well BSD, 

y(x) =Ax + 5tCos(x) +CtCos^(x) +DxCos^(x), 

with Ax=57. 192 kg, 5^=738.415 ke, Cx=2297.880 kg and Dx=557.255 ks. These parameters 
describe a double well with minima at x=180 1 98.8 deg., and a barrier height of about 9 kJ/mol 
at x=180 deg. The barrier height at x=0 deg. is of 30 kJ/mol, out of the range of energies explored 
in these simulations. 

The intermolecular potential is of the non-bonded atom-atom Lennard- Jones type, with pa- 
rameters £=1.70 kJ/mol and o=3.39 A j^Q]. The cut-off radius of the interactions is 15 A and 
correction terms, to the energy and pressure, due to this finite cut-off are taken into account by in- 
tegrating the contribution of an uniform distribution of atoms for distances larger than our cut-off 
of 15 A. 



IV. MOLECULAR DYNAMICS SIMULATIONS 



The MD simulations and samples are entirely identical to those previously performed in refs. 
. The bond constraints are held fixed with the SHAKE algorithm ll22ll and the temperature 



is maintained using the Nose-Hoover chains method ll23ll . although the same behavior was found 
using the standard Nose method I^jQ] in our NPT ensemble simulations. 

Taking the equilibrated samples of a- and Cx'-Sg from our simulations in refs. flQ] as starting 
point, we first did a careful measurement of their crystalline structures (i. e. the lattice parameters 
and volume Vb) by averaging the calculated values over long runs performed in the NPT ensemble. 
The time step of the simulations was of 0.01 ps, the thermalization of the samples was of 80000 
time steps (800 ps) and they were measured in the following 30000 (300 ps) time steps. The lattice 
parameters of each one of the four samples, together with their averaged atomic positions, were 
taken as the equilibrium atomic locations for the corresponding Einstein crystals. 

The next series of MD simulations, as a function of the X parameter (0 < A, < 1), were per- 
formed at constant volume, with a standard time step of 0.005 ps and runs of 160000 time steps 
(800 ps) of thermalization followed by a 40000 time steps (200 ps) of a free trajectory in the phase 
space, that were used to measure the systems. 

For X values close to 1, we had to reduce the time step up to a minimum of 0.0001 ps in 
order to obtain the same accuracy in the total energy as for the other values of X. Accordingly, 
the total number of time steps were increased so as to accumulate nearly the same total time of 
thermalization and storage of the other cases. 

As in refs. 11181 12611 . ten values of the parameter X are defined by a ten point Gauss-Legendre 
quadrature method, used to resolve the thermodynamic integration on the X coordinate (last term 
en eq. 131). 



V. RESULTS 



A. Thermodynamic integration method 

Fig. 121 shows the values of UfreeCX) = {UEinstein — U)-^, obtained with averages over the lasts 
200 ps of each run, for a- and a'-Sg at 300 K. Similar curves Ufreei'X) are calculated for a— and 
a'-5'8 at lOOK, except that for X ~1 < Ufree > is about —260 kJ/mol. For both temperatures the 
values for the a'-Ss sample are systematically lower than those of (X-Sg, the difference increases 
for increasing values of X. 

Note that the values Ufree{X) near A, ~1 are those that weight more in the calculation of the 
integral in eq. IH X ~1 implies strong springs of Einstein crystals and weak interactions of the 



original system, in this case the slope of < UfreeO^ ~ 1) > has a large increment because each 
atom is subjected almost only to the harmonic force of the Einstein crystal, allowing molecular 
configurations very different from those of the real system. Each atom shows a mean square 
displacement that depends only on temperature and not in the forces of the real system. At 
higher temperatures the interactions are more repulsive (because of the shorter mean interatomic 
distances) and < U free > takes more negative values (U has a negative sign in the integrand of 
eq. HI). This can be verified noting that for T = 100 K the minimum < Ufree > is approximately 
-260 kJ/mol, almost half the value at 300 K (-510 kJ/mol, see figure^). 




Figure 2: {Uf,-ee) for a— and a'-S^ at 300K. Each point is a MD run. This curve is the integrand of the 
last term in equation ID from which the free energy is obtained using the Frenkel Ladd thermodynamic 
integration scheme . The lines are a guide to the eyes and were calculated with the Akima interpolation 
algorithm. 



The contributions of the different terms to the total free energy of eq. IH in adimensional units, 
are included in Table HI The first three rows present the same values for all the cases, the first two 
of them are large terms that almost cancel each other. The next three rows show the terms that 
are calculated different in each case, and are obtained from the MD simulation data. The main 
contributions are from the Uq static term and from the thermodynamic integration term. In the last 



row the adimensional free energy values for all the cases are shown. 



Phase/Temp. 


a-58 (r = mK) 


o!-S^{T = imK) 


a-58 {T = 3007*:) 


a'-58 {T = 3007*:) 


3jy_lnA* 


-1162.78 


-1162.78 


-1162.78 


-1162.78 


3(W-l)j , 271 ^ 
2N„„i ^"Vpa*^ 


1148.94 


1148.94 


1148.94 


1148.94 


31n(A') 
2N,„„i 


-0.04032 


-0.04032 


-0.04032 


-0.04032 


Pf/o 


-143.05 


-146.08 


-46.91 


-48.47 




-0.02529 


-0.02529 


-0.02540 


-0.02535 




25.1(4) 


25.5(4) 


21.3(4) 


22.0(4) 




-143.79 


-146.63 


-47.65 


-49.22 


Nmol 


-131.41 


-133.97 


-39.05 


-39.91 



Table I: Contribution of the different terms to the free energy value, calculated by the thermodynamic 
integration method for both temperatures and phases. For the sake of comparison each term has been 
adimensionalized using the definitions A* = A/a, a* = aa^ and Vg* = Vq/ct^ when needed. 



The integral for thermodynamic integration in eq. IHwas also calculated using numerical inte- 
gration for the Akima interpolation curve (see figure|2l) and the difference with the Gauss-Legendre 
quadrature is within 6.4 %. The error of this term, included in tablelHwas obtained with the same 
Gauss-Legendre quadrature using the values UfreeQ^) + Ai7 (X) and UfreeQ^) — (X), where MJ 
is the statistical error of each MD simulation. 

There is a further comment about this free energy calculation method and its use in systems with 
multiple constraints, as are the very commonly used semiflexible molecular models with bond 
length constraints. In our simulations the constraints were not only applied to the real systems 
but also to the corresponding Einstein crystals. In this way we computed free energy differences 
between two constrained systems (the real and the reference) for each one of the four studied 
cases. This is irrelevant as long as we compare the free energy between the a- and a'-Sg crystals 
or between them and their reference systems. However, when computing the absolute free energy 
of each crystal, it has to be taken into account that the bond length constraints affect its MD 
trajectory in the phase space . The correction terms for systems with multiple types of constraints 
were recently discussed by Otter and Briels ll27h . they found that the bond constraint contribution 
is not negligible only when comparing systems of different bond lengths 



B. Quasi-harmonic approximation 



The absolute values of the free energies of the four samples were also calculated via the con- 
tribution of the density of vibrational modes D{v) to the entropy and zero point energy of the 
sample. With the measurement of an unique run in the NPT ensemble, for each one of the four 
real systems, we obtained the values included in Table HTl The samples are equilibrated and mea- 
sured at a pressure of 0. jl0.03 kbar and the pV term is not included to allow a direct comparison 
with the thermodynamic integration method. The volume per molecule configurational en- 
ergy (Uconf), zero point energy (Eq^^), entropy (5) and free energy (F) are shown. The free energy 
calculated with this approximation for various temperatures between 100 and 300 K, not included 
here, showed a systematically lower values for a'-Ss. 



Phase 


a-Ss {T = lOOK) 


a'-S^iT = lOOK) 


a-Ss {T = 300K) 


a'-Ss {T = 2>mK) 


(K,)(A3) 


196.8(3) 


196.0(3) 


203.1(4) 


201.2(2) 


(t/con/) (kJ/mol) 


-114.6(1) 


-116.7(1) 


-105.5(3) 


-108.96(3) 


£e„(kJ/mol) 


16.083 


16.285 


15.950 


15.836 


r5(kJ/mol) 


4.3929 


4.2900 


26.6491 


26.9034 


f(kJ/mol) 


-102.4(3) 


-103.9(3) 


-123.2(1.3) 


-126.0(1.2) 



Table II: Different quantities involved in the free energy calculations via the density of vibrational modes, 
stands for the zero point energy term. 

C. Comparison of both methods 

Table Hill includes the calculated absolute values of free energy, in kJ/mol units, for a— and 
a' — at 100 and 300 K, using both methods. 

We expect higher accuracy in the absolute free energy values for the thermodynamic integration 
method, mainly at higher temperatures when the harmonic approximation is worst than in the case 
of 100 K, due to the fact that anharmonic effects are expected to increase at higher temperatures. 

With the thermodynamic integration method the difference of free energies is calculated about 
^a-a'iT = 300^) ~ 2^ , being a'-Sg the more stable phase. Although the absolute free en- 
ergy values are not similar, the difference between both phases is also of about 2 kJ/mol for the 



T=100K 


F„(kJ/mol) 


F2//(kJ/mol) 


T = 300K 


Fr/(kJ/mol) 


FQff(kJ/mol) 




-109.3(4) 


-102.4(3) 


a-Ss. 


-130.0(1.0) 


-123.2(1.3) 




-111.(3) 


-103.9(3) 




-132.4(1.0) 


-126.0(1.2) 



Table III: Comparison of the absolute free energy values using thermodynamic integration (TI) and Quasi- 
harmonic approximation (QH). 

quasi-harmonic method. This difference can be experimentally measured and closely follows that 
between the corresponding Uq and Uconf included in table II VI 



Phase 


T{K) 


Fr/(kJ/mol) 


Fe//(kJ/mol) 


UoikJ/mol) 


Uconf{kJ /mol) 


a-Ss 


100.3(1.9) 


-109.3(4) 


-102.4(3) 


-118.936 


-114.6(1) 


a'-Ss 


100.8(1.9) 


-111.39(34) 


-103.94(25) 


-121.461 


-116.7(1) 


a-Ss 


299.6(5.1) 


-130.27(96) 


-123.2(1.3) 


-117.008 


-105.5(3) 




300.9(4.9) 


-132.4(1.0) 


-126.0(1.2) 


-120.901 


-108.0(3) 



Table IV: Comparison of the free energy values for the thermodynamic integration and the quasi-harmonic 
approximation methods. The values of the energy for the mean positions {Uo) and the mean configuration 
energy (Uconf) are also shown for the two phases and the two studied temperatures. 



VI. CONCLUSIONS 

In this paper we calculated the free energy of orthorhombic a-Sg crystals and compared it with 

nnnn 

that of monoclinic a'-Sg laLZllilllU], a previously proposed phase that is calculated with lower 
mean volume per molecule and lower potential energy than a-Sg crystals. The monoclinic a'- 
Sg crystal was obtained in the course of MD simulations, using simple inter- and intramolecular 
potential models, but has not been experimentally observed Therefore, in this work we 

performed an accurate measurement of the free energy difference between both crystals at 100 
and 300 K (above and below the spontaneous transition found in refj^), as given by the simple 
molecular model of |Q]. The thermodynamic integration method fl, IJ] is implemented, taking as 
reference system the corresponding Einstein crystals. For the sake of comparison we also included 
calculations of the free energy in the quasi-harmonic approximation, a fast estimation method that 
can be easily implemented in MD simulations. 



The accurate free energy calculations performed here still show that monoclinic a'-Sg is calcu- 
£than the orthorhombic experimentally observed (X-Sg, in accord with our previous 
performed with the same simple molecular model and in the quasi-harmonic ap- 



lated more stab 



calculations 
proximation. 

As regards the accuracy of the free energy calculations, we found that for the thermodynamic 
integration method the results are extremely sensitive to the accurate calculation of the values 
U (A,) . In particular for ?i ~ 1 long MD simulations are needed. Although the absolute free energy 
values differ for both calculation methods, both measure a difference of about 2 kJ/mol between 
a- and a'-Sg crystals, being a'-Sg the more stable. 

The approximated method for calculating the free energy via the contribution of the vibrational 
modes to the energy and entropy of the system turns out to be very useful to determine relative 
values between polymorphic phases or between different temperatures of a given sample. This 
result is very important, taking into account that this method is, at least, ten times faster than the 
thermodynamic integration method, that requires at least ten lengthly MD simulations (one MD 
run for each value of the X parameter), in order to calculate the free energy of one sample at a 
given temperature. 
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I. APPENDIX 

We estimate here the difference in the calculated statistical averages for MD simulations per- 
formed with and without molecular bond length constraints. This topic has been extensively 
treated in literature and is not of minor importance Indeed, it has been widely 

studied in the case of free energy calculations that involve a reaction coordinate 

We recall the results of chapter 10 from ref. [16i | and the appendix from ref. 112 811 . The ex- 
pression for the ratio of molecular position probability for constrained and unconstrained system 
is: 



p(qi,...,q3iv) 



(1) 

P(qi,-,q3yv-/,qi =Oi,-,q/ =<7/J 
where we explicitly note the constrained degrees of freedom with q^' = O/, being the standard 

bond length constraint, 



J = 0), 



(2) 



and / constrained bonds are assumed. The I x I Z matrix, in the right hand side of[I]is defined 



as M 



281 



-'mn 



(3) 



and 1 1 stands for the determinant. The origin of Z is from the fact that not only air;) = must be 
required to constraint the bonds in the simulations, but also o(r,) = at all times OOll . It should 
be noted that q in eq. [T]is used for generalized coordinates while for cartesian ones. 

For the Sg molecule we have 24 degrees of freedom, from which 8 have been constrained in our 
model (see sec. HiHi . Therefore we have 16 degrees of freedom with 8 bond constraints following 
eq. El with d =2.06 A. Matrix Z can be explicitly calculated for a molecule from eq. 13 

(pl8)^ 



eq. El with d =2.06 A. Matrix 

/ 2 cos(|3i2) 
cos((3i2) 2 C0S((323 
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C0S((334) 
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C0S((345) 
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C0S(|356) 











C0S((356) 

2 

COs(p67) 











cos(p67; 
2 

C0S(|378) 



cos( 






C0S((378) 

2 



get: 



where |3,j is the bending bond angle between consecutive angles z, j and j ' + 1. 
If we take a mean (3 value for all the bending angles and evaluate the determinant in eq. d we 



|Z| =256-576cos(|3)2 + 400cos(p)4-88cos(|3)^ (4) 

We can now take a statistical mean value over the NVT ensemble, averaging eq. (Hover a 
MD run and normalizing with the maximum value. This provides an estimation of the effect 



of constraints over statistical averages such as used in eq. IH to carry out the thermodynamic 
integration. We got the following main difference between the constrained and unconstrained 
probabilities atT = 300K: 

' -0.954. 



Pconstrained 

The closeness to 1 indicates that the effect of bond constrains is practically negligible for sim- 
ulations in the studied range of thermodynamic variables. From this estimation we can conclude 
that the points obtained in figure El are almost identical to those that can be obtained in a similar 
simulation without bond constrains, that would have been quite more expensive in CPU time. This 
simple calculation, that we present here for close ring-molecules, is complementary to that already 
obtained in ref. 113 ih for open chain molecules. 



Contribution of stretching modes to free energy 

We include here an estimation of the free energy difference between samples with and without 
molecular bond constrains, under the framework of quasi-harmonic approximation commented in 
sec. IllBI If we think in the normal modes as harmonic oscillators of frequency v,, we have a 
'trivial' partition function for each oscillator of the form lni : 



1 - e-'''^' 



e 

5=0 



where h is the Planck constant, k the Boltzmann factor and s the positive integer quantum 
number of the oscillator of frequency v,. For a system of such independent oscillators we have a 
partition function: 

Qsunul = l Z)(v)(^--p^)JV, (5) 

where D{\) is the density of vibrational states already defined in sec. Ill Bh nd Qswmi is cal- 
culated from our bond constrained MD simulation. The partition function of the unconstrained 
system Q can be estimated in the following way: taking into account that for the Sg molecule the 
stretching vibrational modes are well above bending and torsional modes in the v scale (see sec. 
imT) . we have separable contributions for the total partition function of the molecule and 

Q = / D,i^ul{y){- R7-) + / D,tretchMi- oT-) dV = Qsimul + Qstretch, (6) 

Jo l-e P'^^ Jv„,flv 1 - 



where we have assumed that the interval {0^v,„ax) accounts for the bending, torsional (dihedral) 
and lattice vibration modes and the high frequencies of bond stretching are 'isolated' in the interval 
(Vmax, °°) ■ The first term in eq. |6lis obtained from our simulation and the second one can be simply 
approximated by: 

Qstretch = N'"''''^- 



1 _ e-^hV„re,ch ' ' 

here M'^"'^"^'^ are the 8 stretching degrees of freedom (constrained in the MD simulations) and 
^stretch is the mean stretching frequency. This is equivalent to use the relationship Dstretch{^) ~ 
^stretch ^ g^-y _ y^f^gj^/,). The free energy difference between considering or not the stretching 
degrees of freedom can be estimated, under the quasi-harmonic approximation by: 

PAF = mot-Fsimui) = ln( ^7"^ ) 

\lsimul "r stretch 

For T = 300K and the sulphur mean stretching frequency \ stretch = 453. 8 cm^^ we get |3AF = 
0.075 which is a relatively and also an absolute small number when compared with the values 
of table D 

This simple calculation can be considered complementary to those of ref. 
provided a thorough discussion of the role of bond constrains in free energy calculations. 
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